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LONG-TERM  GOALS 

The  long-tenn  technical  goal  is  to  design,  develop  and  test  the  next  generation,  primitive  equation 
ocean  model  for  high-resolution  scientific  (ROMS:  Regional  Ocean  Modeling  System)  and  operational 
(TOMS:  Terrain- following  Ocean  Modeling  System)  applications.  Currently,  both  modeling  systems 
are  identical  and  extensively  used  by  the  research  and  operational  communities  around  the  world.  Our 
aim  is  to  produce  an  open-source,  terrain-following,  ocean  community  model  for  regional  nowcasting 
and  forecasting  that  includes  advanced  data  assimilation,  ensemble  prediction,  and  analysis  tools  for 
adaptive  sampling  and  circulation  dynamics/stability/sensitivity,  which  is  highly  relevant  to  ONR 
objectives.  This  project  will  improve  the  ocean  modeling  capabilities  of  the  U.S.  Navy  for  relocatable, 
coastal,  coupled  atmosphere-ocean  forecasting  applications. 

OBJECTIVES 

The  main  objectives  of  this  project  are: 

•  To  develop  and  test  a  robust  ocean  modeling  framework  for  relocatable  coastal  ocean 
prediction  applications. 

•  To  develop  algorithms  and  tools  for  improving  U.S.  Navy  coupled  atmosphere-ocean 
forecasting  capabilities. 

•  To  develop  advanced  4-dimension  Variational  (4D-Var)  data  assimilation  capabilities  and 
analysis  algorithms  for  observation  sensitivity,  observation  impact,  adaptive  sampling,  and 
forecast  errors  and  uncertainties. 

•  To  develop  adjoint-based  ocean  prediction  analysis  tools  similar  to  those  available  in 
Numerical  Weather  Prediction  (NWP)  for  the  atmosphere  for  circulation  stability,  sensitivity 
analysis,  and  ensemble  prediction. 

•  To  develop  multiple  grid  nesting  capabilities  to  resolve  complex  geographical  regions  and 
circulation  regimes. 
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•  To  develop  the  framework  for  multiple  model  (atmosphere,  sea-ice,  waves)  coupling  using 
available  libraries. 

•  To  provide  the  ocean  modeling  community  with  the  current  state-of-the-art  knowledge  in 
dynamics,  numerical  schemes,  and  computational  algorithms  technology.  ROMS  is  freely 
distributed  (www.myroms.org)  to  the  Earth’s  modeling  community  and  has  thousands  of  users 
worldwide. 

•  To  engage  the  ocean  modeling  community  by  organizing  annual  scientific  workshops  and 
training. 

APPROACH 

ROMS/TOMS  is  a  three-dimensional,  free-surface,  terrain-following  ocean  model  that  solves  the 
Reynolds-averaged  Navier-Stokes  equations  using  the  hydrostatic  vertical  momentum  balance  and 
Boussinesq  approximation  (Haidvogel  et  al.  2000,  2008;  Shchepetkin  and  McWilliams,  2005,  2009). 
The  governing  dynamical  equations  are  discretized  on  a  vertical  coordinate  that  depends  on  the  local 
water  depth.  The  horizontal  coordinates  are  orthogonal  and  curvilinear  allowing  Cartesian,  spherical, 
and  polar  spatial  discretization  on  an  Arakawa  C-grid.  Its  dynamical  kernel  includes  accurate  and 
efficient  algorithms  for  time-stepping,  advection,  pressure  gradient  (Shchepetkin  and  McWilliams 
2003,  2005),  several  subgridscale  parameterizations  (Durski  et  al.,  2004;  Warner  et  al.,  2005)  to 
represent  small-scale  turbulent  processes  at  the  dissipation  level,  and  various  bottom  boundary  layer 
formulations  to  determine  the  stress  exerted  on  the  flow  by  the  bottom. 

Several  adjoint-based  algorithms  exist  to  explore  the  factors  that  limit  the  predictability  of  the 
circulation  in  regional  applications  for  a  variety  of  dynamical  regimes  (Moore  et  al.,  2004,  2009). 

These  algorithms  use  the  ideas  of  Generalized  Stability  Theory  (GST)  in  order  to  identify  the  most 
unstable  directions  of  state-space  in  which  errors  and  uncertainties  are  likely  to  grow.  The  resulting 
singular  vectors  can  be  used  to  construct  ensembles  of  forecasts  by  perturbing  initial  and  boundary 
conditions  (optimal  perturbations)  and/or  surface  forcing  (stochastic  optimals).  Perturbing  the  system 
along  the  most  unstable  directions  to  the  state-space  yields  information  about  the  first  (ensemble  mean) 
and  second  (ensemble  spread)  moments  of  the  probability  density  function.  Given  an  appropriate 
forecast  skill  measure,  the  circulation  is  predictable  if  low  spread  and  unpredictable  if  large  spread. 

ROMS/TOMS  uniquely  supports  three  different  4D-Var  data  assimilation  methodologies  (Moore  et  al., 
2011a,  b):  a  primal  fonn  of  the  incremental  strong  constraint  4D-Var  (I4D-Var),  a  strong/weak 
constraint  dual  form  of  4D-Var  based  on  the  Physical-space  Statistical  Analysis  System  (4D-PSAS), 
and  a  strong/weak  constraint  dual  form  of  4D-Var  based  on  the  indirect  representer  method  (R4D-Var). 
In  the  dual  formulations,  the  search  for  the  best  ocean  circulation  estimate  is  in  the  subspace  spanned 
only  by  the  observations,  as  opposed  to  the  full  space  spanned  by  the  model  as  in  the  primal 
formulation.  Although  the  primal  and  dual  formulations  yield  identical  estimates  of  the  ocean 
circulation  for  the  same  a  priori  assumptions,  there  are  practical  advantages  and  disadvantages  to  both 
approaches  (Moore  et  al.,  2011a,  b,  c).  To  our  knowledge,  ROMS/TOMS  is  the  only  open-source, 
ocean  community-modeling  framework  supporting  all  these  variational  data  assimilation  methods  and 
other  sophisticated  adjoint-based  algorithms. 
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There  are  several  biogeochemical  models  available  in  ROMS.  In  order  of  increasing  ecological 
complexity  these  include  three  NPZD-type  models  (Franks  et  al.,  1986;  Powell  et  al.,  2006;  Fiechter  et 
al.,  2009),  a  nitrogen-based  ecosystem  model  (Fennel  et  al.,  2006,  2008),  a  Nemuro-type  lower  level 
ecosystem  model  (Kishi  et  al.,  2007),  and  a  bio-optical  model  (Bissett  et  al.,  1999). 

ROMS  includes  a  sediment-transport  model  with  an  unlimited  number  of  user-defined  cohesive  (mud) 
and  non-cohesive  (sand)  sediment  classes  (Warner  et  al,  2008).  Each  class  has  attributes  of  grain 
diameter,  density,  settling  velocity,  critical  stress  threshold  for  erosion,  and  erodibility  constant.  A 
multi-level  bed  framework  tracks  the  distribution  of  every  size  class  in  each  layer  and  stores  bulk 
properties  including  layer  thickness,  porosity,  and  mass,  allowing  the  computation  of  bed  morphology 
and  stratigraphy.  Also  tracked  are  bed-surface  properties  like  active-layer  thickness,  ripple  geometry, 
and  bed  roughness.  Bedload  transport  is  calculated  for  mobile  sediment  classes  in  the  top  layer. 

ROMS  is  a  very  modern  and  modular  code  written  on  F90/F95.  It  uses  C-preprocessing  to  activate  the 
various  physical  and  numerical  options.  The  parallel  framework  is  coarse-grained  with  both  shared- 
memory  (OpenMP)  and  distributed-memory  (MPI)  paradigms  coexisting  in  the  same  code.  Because  of 
its  construction,  the  parallelization  of  the  adjoint  is  only  available  for  MPI.  Several  coding  standards 
have  been  established  to  facilitate  model  readability,  maintenance,  and  portability.  All  the  state  model 
variables  are  dynamically  allocated  and  passed  as  arguments  to  the  computational  routines  via 
dereferenced  pointer  structures.  All  private  arrays  are  automatic;  their  size  is  determined  when  the 
procedure  is  entered.  This  code  structure  facilitates  computations  over  nested  grids  (Warner  et  al., 
2011). 

WORK  COMPLETED 

Three  types  of  nesting  capabilities  have  been  designed  and  coded  in  ROMS:  (i)  refinement  grids  which 
provide  increased  resolution  (3:1,  5:1,  or  7:1)  in  a  specific  region;  (ii)  mosaics  which  connect  several 
grids  along  their  edges,  and  (iii)  composite  grids  which  allow  overlap  regions  of  aligned  and  non- 
aligned  grids.  The  mosaic  and  composite  grid  code  infrastructures  are  identical.  The  differences  are 
geometrical  and  primarily  based  on  the  alignment  between  adjacent  grids.  All  the  mosaic  grids  are 
exactly  aligned  with  the  adjacent  grid.  In  general,  the  mosaic  grids  are  a  special  case  of  the  composite 
grids. 

The  nesting  development  in  ROMS  was  divided  into  three  phases  due  to  its  complexity.  Phase  I 
included  substantial  modifications  of  the  numerical  kernels  to  allow  a  generic  treatment  of  the  spatial 
horizontal  operators  in  the  nesting  contact  regions.  Phase  II  included  an  overhaul  of  ROMS  lateral 
boundary  conditions  to  facilitate,  in  a  generic  way,  their  processing  or  not  in  applications  with  nested 
grids.  Phase  III  included  the  data  managing  and  time-stepping  infrastructure  for  one  or  more  nesting 
layers.  These  phases  are  described  in  detail  in  last  year’s  (FY1 1)  annual  report.  Phase  I  was  released 
to  the  community  as  ROMS  3.5  on  April  25,  2011  whereas  Phase  II  was  released  as  ROMS  3.6  on 
September  23,  2011.  The  coding  of  Phase  III  was  completed  during  this  fiscal  year  cycle  and  is 
currently  under  extensive  testing  before  it  is  released  to  the  user  community. 
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RESULTS 


During  this  fiscal  year  (FY12),  the  primary  focus  was  the  development  and  coding  of  ROMS  nesting 
Phase  III.  Several  modules  were  added  to  process  the  information  that  is  required  in  the  contact 
regions,  what  information  needs  to  be  exchanged  from/to  another  grid,  and  when  to  exchange  it. 

A  contact  region  is  an  extended  section  of  the  grid  that  overlays  or  is  adjacent  to  a  nested  grid.  A 
contact  point  is  a  grid  cell  inside  the  contact  region.  Each  contact  region  has  a  receiver  grid  and  donor 
grid.  The  contact  points  of  the  receiver  grid  are  processed  using  the  donor  grid  cell  containing  the 
contact  point.  The  contact  points  are  coincident  if  the  receiver  and  donor  grids  occupy  the  same 
position.  If  coincident  grids,  the  receiver  grid  contact  point  data  is  just  filled  using  the  donor  data. 
Otherwise  (if  not  coincident)  the  receiver  grid  data  is  linearly  interpolated  from  the  donor  grid  cell 
containing  the  contact  point.  There  is  a  duality  in  ROMS  nested  grids:  data  donor  in  one  contact  region 
and  data  receiver  in  its  conjugate  contact  region.  The  exchange  of  information  is  always  two-way.  This 
explains  why  we  prefer  to  use  the  donor  and  receiver  categories  instead  of  the  usual  parent  and  child 
descriptions  found  in  the  literature. 

The  contact  points  are  processed  outside  of  ROMS  and  all  the  connectivity  within  nested  grids  is  read 
from  an  input  NetCDF  file.  This  facilitates  the  configuration  of  various  grid  Sub-Classes.  It 
tremendously  simplifies  the  processing  of  such  points  in  parallel  computations.  It  also  gives  the  user 
full  editing  control  of  the  contact  points.  The  contact  points  need  to  be  set-up  once  for  a  particular 
application.  Several  Matlab  utilities  were  developed  to  process  ROMS  nested  grids  capabilities.  These 
scripts  are  too  technical  to  be  described  here,  but  detailed  information  can  be  found  in  WikiROMS 
(https://www.myroms.org/wiki/index.php/Matlab_Scripts). 

In  grid  refinement  applications,  the  information  is  exchanged  at  the  beginning  of  the  full  time-step  (top 
of  main2d  or  main3d).  Contrarily,  in  composite  grid  applications,  the  information  is  exchanged 
between  each  sub-time  step  call  in  main2d  or  main3d.  That  is,  the  composite  donor  and  receiver  grid 
need  to  sub-time  step  the  2D  momentum  equations  before  any  of  them  start  solving  and  coupling  the 
3D  momentum  equations.  Therefore,  the  governing  equations  are  solved  and  nested  in  a  synchronous 
fashion.  The  concept  of  nesting  layers  is  introduced  to  allow  applications  with  both  composite  and 
refinement  grid  combinations,  as  shown  in  Figure  3a  for  the  Refinement  and  Partial  Boundary 
Composite  Sub-Class.  Here,  two  nesting  layers  are  required.  In  the  first  nesting  layer,  the  composite 
grids  1  and  2  are  synchronously  sub-time  stepped  for  each  governing  equation  term.  Then,  in  the 
second  and  last  nesting  layer,  the  refinement  grid  3  is  time  stepped  and  the  information  between  donor 
(1)  and  receiver  refinement  (3)  grids  are  exchanged  at  the  beginning  of  the  time  step  for  both  contact 
regions:  coarse-to-fine  and  fine-to-coarse. 

The  ROMS  nested  grid  design  includes  three  Super-Classes  and  several  Sub-Classes: 

1 .  Composite  Grids  Super-Class: 

a.  Mosaic  Grids  Sub-Class 

b.  Composite  Overlap  Grids  Sub-Class 

c.  Complex  Estuary  Composite  Grids  Sub-Class 

d.  Partial  Boundary  Composite  Grids  Sub-Class 

2.  Refinement  Grids  Super-Class: 
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a.  Single  Refinement  Sub-Class 

b.  Multiple  Refinement  Sub-Class 

3.  Composite  and  Refinement  Combination  Super-Class: 

a.  Refinement  and  Partial  Boundary  Composite  Sub-Class 

b.  Complex  Estuary  Refinement-Composite  Sub-Class 

Hence,  there  are  several  possibilities  and  combinations.  The  design  is  flexible  enough  to  allow 
complex  nested  grid  configurations  in  coastal  applications.  These  capabilities  are  better  illustrated  with 
diagrams  of  the  different  Sub-Classes  shown  in  Figures  1,  2,  and  3.  The  diagrams  are  all  in  terms  of 
computational  coordinates  (£,,  q)  or  fractional  (/,/)  coordinates.  Therefore  we  have  squares  and 
rectangles  that  may  map  to  physical  curvilinear  coordinates.  Here,  Ngrids  is  the  number  of  nested 
grids,  NestLayers  is  the  number  of  nested  layers,  GridsInLayer  is  a  vector  with  the  number  of  grids 
in  each  nested  layer,  and  Ncontact  is  the  number  of  nesting  contact  regions,  (Ngrids-1)*2. 
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a) 


c) 


Mosaic  Sub-Class: 

Ngrids  =  2 
NestedLayers  =  1 
GridsInLayer  =  2 
Ncontact  =  2 


b) 

1 

2 

Composite  Overlap 
Sub-Class: 

Ngrids  =  2 


Complex  Estuary 
Composite  Sub-Class: 

Ngrids  =  6 
NestedLayers  =  1 
GridsInLayer  =  6 
Ncontact  =10 


Partial  Boundary 
Composite  Sub-Class: 

Ngrids  =  2 
NestedLayers  =  1 
GridsInLayer  =  2 
Ncontact  =  2 


Figure  1:  Composite  Grids  Super-Class:  a)  Mosaic  Sub-Class,  b)  Composite  Overlap  Sub-Class,  c) 
Complex  Estuary  Composite  Sub-Class,  and  d)  Partial  Boundary  Composite  Sub-Class.  Each 
diagram  shows  the  number  of  nested  grids  (Ngrids),  number  of  nested  layers  (NestedLayers), 
number  of  grids  in  each  nested  layer  (GridsInLayer),  and  number  of  contact  regions  (Ncontact). 
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b) 


Single  Refinement 
Sub-Class: 

Ngrids  =  2 
NestedLayers  =  2 
GridsInLayer  =  1  1 
N  contact  =  2 


Multiple  Refinement 
Sub-Class: 

Ngrids  =  4 
NestedLayers  =  3 
GridsInLayer  =12  1 
Ncontact  =  6 


Figure  2:  Refinement  Super-Class:  a)  Single  Refinement  Sub-Class  and  b)  Multiple  Refinement 
Sub-Class.  Each  diagram  shows  the  number  of  nested  grids  (Ngrids),  number  of  nested  layers 
(NestedLayers),  number  of  grids  in  each  nested  layer  (GridsInLayer),  and  number  of 

contact  regions  (Ncontact). 
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Refinement  and  Partial 
Boundary  Composite 
Sub-Class: 

Ngrids  =  3 
Nested  Layers  =  2 
GridsInLayer  =  2  1 
Ncontact  =  4 


Complex  Estuary 

Refinement-Composite 

Sub-Class: 

Ngrids  =  3 
Nested  Layers  =  2 
GridsInLayer  =  1  2 
Ncontact  =  4 


Figure  3:  Composite  and  Refinement  Combination  Super-Class:  a)  Refinement  and  Partial 
Boundary  Composite  Sub-Class  and  b)  Complex  Estuary  Refinement-Composite  Sub-Class.  Each 
diagram  shows  the  number  of  nested  grids  (Ngrids),  number  of  nested  layers  (NestedLayers), 
number  of  grids  in  each  nested  layer  (GridsInLayer),  and  number  of  contact  regions  (Ncontact). 

Figure  4  shows  a  realistic  nesting  configuration  in  a  US  East  Coast  application  using  the  Complex 
Estuary  Refinement-Composite  Sub-Class.  It  illustrates  ROMS  nesting  capabilities  with  complex 
coastlines  and  estuaries.  The  coarser  grid,  ESPRESSO  (130  x  82),  has  an  average  resolution  of 
dx=7.5km,  dy=5.8km.  The  nested  grids  (p-points  mesh)  are  color  coded  for  convenience  to  show  the 
strategy  used  to  better  resolve  the  Delaware  and  Chesapeake  Estuary  Systems.  The  red  and  green  are 
refinement  grids  whereas  blue  and  magenta  are  composite  grids.  The  refinement  ratio  is  1:7.  An 
intermediate  1 :7  refinement  grid  is  created  using  Matlab  script  coarse2fine.m  that  included  both  the 
Delaware  and  Chesapeake  Estuary  Systems.  Then,  the  Matlab  script  gridextract.m  is  used  to  extract 
the  Delaware  Bay  refinement  grid  (58  x  142)  and  Delaware  River  composite  grid  (42  x  55).  Similarly, 
grid  extr act. m  is  used  to  extract  the  Chesapeake  Bay  outer  refinement  grid  (135  x  142)  and 
Chesapeake  Bay  inner  composite  grid  (233  x  212). 
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Figure  4:  A  US  east  coast  example  of  the  Complex  Estuary  Refinement-Composite  Sub-Class:  a) 
full  nested  grid  configuration,  b)  zoom  of  the  Delaware  and  Chesapeake  estuaries  refinement  (green 
and  red)  and  composite  (magenta  and  blue)  grids,  c)  further  zoom  of  Delaware  estuary  contact 
region  between  refinement  and  composite  connected  grids,  and  d)  further  zoom  of  Chesapeake 
estuary  contact  region  between  refinement  and  composite  connected  grids. 


9 


The  next  user  workshop  will  be  held  in  Rio  de  Janeiro,  Brazil,  October  22-24,  2012.  A  special  fourth 
day  (October  25)  has  been  added  to  focus  modern  observational,  modeling,  and  data  assimilation 
systems,  which  includes  special  lectures  and  discussions.  Over  100  participants  are  expected  this  year; 
far  more  than  at  any  previous  workshop. 

IMPACT/APPLICATIONS 

This  project  will  provide  the  ocean  modeling  community  with  a  freely  accessible,  well  documented, 
open-source,  terrain-following,  ocean  model  for  regional  nowcasting  and  forecasting  that  includes 
advanced  data  assimilation,  ensemble  prediction,  and  analysis  tools  for  adaptive  sampling  and 
circulation  dynamics,  stability,  and  sensitivity. 

TRANSITIONS 

The  full  transition  of  ROMS/TOMS  to  the  operational  community  is  likely  to  occur  in  the  future. 
However,  the  ROMS/TOMS  algorithms  are  now  available  to  the  developers  and  scientific  and 
operational  communities  through  the  website  http://www.myroms.org/. 

RELATED  PROJECTS 

The  work  reported  here  is  related  to  other  already  funded  ONR  projects  using  ROMS.  In  particular,  the 
PI  (H.  Arango)  closely  collaborates  with  A.  Moore  (adjoint-based  algorithms)  at  University  of 
California,  Santa  Cruz,  A.  Miller  and  B.  Comuelle  (ROMS  adjoint  and  variational  data  assimilation)  at 
Scripps  Institute  of  Oceanography,  E.  Di  Lorenzo  (Southern  California  predictability)  at  Georgia 
Institute  of  Oceanography,  and  J.  Wilkin  (Mid-Atlantic  Bight  variational  data  assimilation)  at  Rutgers 
University. 
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